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ABSTRACT 

Recent Chandra and XMM-Newton observations of a number of X-ray "dim" pul- 
sating neutron stars revealed quite unexpected features in the emission from these 
sources. Their soft thermal spectrum, believed to originate directly from the star sur- 
face, shows evidence for a phase- varying absorption line at some hundred eVs. The 
pulse modulation is relatively large (pulsed fractions in the range ~ 12%-35%), the 
pulse shape is often non-sinusoidal, and the hard X-ray color appears to be anti- 
correlated in phase with the total emission. Moreover, the prototype of this class, RX 
J0720. 4-3125, has been found to undergo rather sensible changes both in its spectral 
and timing properties over a timescale of a few years. All these new findings seem dif- 
ficult to reconcile with the standard picture of a cooling neutron star endowed with a 
purely dipolar magnetic field, at least if surface emission is produced in an atmosphere 
on top of the crust. In this paper we explore how a dipolar -|-quadrupolar star-centered 
field influence the properties of the observed lightcurves. The phase-resolved spectrum 
has been evaluated accounting for both radiative transfer in a magnetized atmosphere 
and general relativistic ray-bending. We computed over 78000 lightcurves varying the 
quadrupolar components and the viewing geometry. A comparison of the data with 
our model indicate that higher order multipoles are required to reproduce the obser- 
vations. 

Key words: Radiative transfer — stars: neutron — pulsars: general — X-rays: stars 
— stars: individual (RX J0720.4-3125, RX J0420.0-5022, RX J0806.4-4123, RBS 1223, 
RBS 1774) 



1 INTRODUCTION 

Over the last few years a number of high resolution spectral 
and timing observations of thermally emitting neutron stars 
(NSs) have become available thanks to new generation X-ray 
satellites (both Chandra and XMM-Newton) , opening new 
perspectives in the study of these sources. Thermal emis- 
sion from isolated NSs is presently observed in more than 
20 sources, including active radio pulsars, soft 7-repeaters, 
anomalous X-ray pulsars, Geminga and Geminga-like ob- 
jects, and X-ray dim radio-silent NSs. There is by now a 
wide consensus that the soft, thermal component directly 
originates from the surface layers as the star cools down. 
If properly exploited, the information it conveys are bound 
to reveal much about the physics of neutron stars, shedding 
light on their thermal and magnetic surface distribution and 
ultimately probing the equation of state of matter at supra- 
nuclear densities. 



Although thermal surface emission seems indeed to be 
an ubiquitous feature in isolated NSs, a power-law, non- 
thermal component (likely produced in the star magneto- 
sphere) is present in most sources, where it often dominates 
the X-ray spectrum. Moreover, the intrinsic X-ray emission 
from young radio-pulsars may be significantly contaminated 
by the contribution of the surrounding supernova remnant. 
In this respect the seven dim X-ray sources discovered by 
ROSAT (hereafter XDINSs) are a most notable exception. 
In a sense, one may claim that these are the only "genuinely 
isolated" NSs and their soft thermal emission is unmarred 
by (non-thermal) magnetospheric activity nor by the pres- 
ence of a sup ernova rem nant or a bin ary compan ion (see e.g. 
iTreves et allEoOa and .HaberJl2004l for reviews; IZane et alJ 
|2005f) . XDINSs play a key role in compact objects astro- 
physics: these are the only sources in which we can have a 
clean view of the compact star surface, and as such offer an 
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unprecedented opportunity to confront theoretical models 
of neutron star surface emission with observations. 

The XDINSs X-ray spectrum is with no exception 
blackbody-like with temperatures in the range ~ 40-100 eV 
and, thus far, pulsations have been detected in five sources, 
with periods in the range 3-11 s (see Table and refs. 
therein) . In each of the five cases the pulsed fraction is rela- 
tively large (--^ 12%-35%). Quite surprisingly, and contrary 
to what one would expect in a simple dipolar geometry, of- 
ten the hardness rat i o is minimum at t he pulse maximum 
jCropper et alj|200lt iHaberl et allbOOaH . Broad absorption 
features have been detected around ~ 300-700 eV in all pul- 
sating XDINSs and the line strength appears to vary with 
the pulse phase. In addition, the X-ray light curves exhibit a 
certain asymmetry, with marked deviations fr om a pure si- 
nusoidal shape at least in the case of RBS 1223 iHaberl et alJ 
L2003; SchwoDC ct al. 2005). XDINSs were unanimously be- 
lieved to be steady sources, as indicated by several years of 
observations for the brightest of them. Unexpectedly, and 
for the first ever time, XMM-Newton observations have re- 
cently revealed a substantial change in the spectral shape 
and pulse profile of the second most lumin ous source, RX 
J0720. 4-3125, over a timescale of ~ 2 yr llDe Vries et alJ 
l2004l:IVink et alJl2004^ . Possible variations in the pulse pro- 
file of RX J0420.0-5022 over a similar timescale (~ 0.5 yr) 
have also been reported, although only at a low significance 
level (iHaberl et al.ll2004l) . 

In the standard picture, emission from an isolated, cool- 
ing NS arises when thermal radiation originating in the out- 
ermost surface layers traverses the atmosphere which covers 
the star crust. Although the emerging spectrum is thermal, 
it is not a blackbody because of radiative transfer in the 
magnetized atmosphere and the inhomogeneous surface tem- 
perature distribution. The latter is controlled by the crustal 
magnetic field, since thermal conductivity across the field is 
highly suppressed, and bears the imprint of the field topol- 
ogy. Besides the spectrum, radiative transfer and the sur- 
face temperature distribution act together in shaping the 
X-ray lightcurve. Pulse profiles produced by the thermal 
surface distribution induced by a simple core-centered dipo- 
lar ma gnetic field have been investigated long ago bv lPage| 
il995() . under the assumption that each surface patch emits 
(isotropic) blackbody radiation. Because of gravitational ef- 
fects and of the smooth temperature distribution (the tem- 
perature monotonically decreases from the poles to the equa- 
tor), the pulse modulation is quite modest (pulsed frac- 
tion < 10%) for reasonable values of the star radius. More- 
over, being the temperature distribution symmetric about 
the magnetic equator, the pulse shape itself is always sym- 
metrical, regardless of the viewing geometry. Larger pulsed 
fractions may be reached by the proper inclusion of an at- 
mosphere. In fact, in a strongly magnetized medium photon 
propagation is anysotropic and o ccurs preferentially along 
the field (magnetic beaming, e.g. IPavlov et aljIlQM) . Nev- 
ertheless, retaining a dipolar temperature distribution will 
always result in a symmetric pulse profile. 

The quite large pulsed fraction, pulse asymmetry, and 
possibly long-term variations, recently observed in XDINSs 
seem therefore difficult to explain by assuming that the ther- 
mal emission originates at the NS surface, at least when 
assuming that the thermal surface distribution is that in- 
duced by a simple core-centered dipolar magnetic field. It 



should be stressed that, although the dipole field is a con- 
venient approximation, the real structure of NSs magnetic 
field is far from been understood, e.g. it is still unclear if the 
field threads the entire star or is co nfined in the crust only 
fe.g. iGeppert. Kiiker fc Pagdl2004l and references therein). 
Whatever the case, there are both observational and theo- 
retical indications that the NS surface field is "patch y" (e.g. 
iGeppert. Rheinhardt fc G if 2003': lUrpin fc Gil2004l and ref- 
erences therein). The effects of a more complex field g eome- 
try have been investigated bv lPage fc SarmientoHl99(jl . who 
considered a star-centered dipole-|-quadrupole field, again 
assuming isotropic emission. The presence of multipolar 
components induces large temperature variations even be- 
tween nearby regions and this results in larger pulsed frac- 
tions and asymmetric pulse profiles. 

The high quality data now available for thermally emit- 
ting NSs, and XDINSs in particular, demand for a detailed 
modelling of surface emission to be exploited to a full ex- 
tent. Such a treatment should combine both an accurate 
formulation of radiation transport in the magnetized atmo- 
sphere and a quite general description of the thermal and 
magnetic surface distributions, which, necessary, must go 
beyond the simple dipole approximation. The ultimate goal 
is to produce a completely self-consistent model, capable to 
reproduce simultaneously both the spectral and timing prop- 
erties. In this paper we take a first step in this direction 
and present a systematic study of X-ray lightcurves from 
cooling NSs, accounting for both a quadrupolar magnetic 
field (in addition to the core-centered dipole) and radia- 
tive transfer in the magnetized atmosphere. We computed 
over 78000 model lightcurves, exploring the entire parame- 
ter space, both in the geometrical angles and the quadrupo- 
lar components. This large dateiset has been analyzed using 
multivariate statistical methods (the principal component 
analysis and the cluster analysis) and we show that a non- 
vanishing quadrupolar field is required to reproduce the ob- 
served XDINS pulse profiles. 



2 THE MODEL 

2.1 Going Quadrupolar 

In this section we describe the approach we use to compute 
the phase-dependent spectrum emitted by a cooling NS as 
seen by a distant observer. This issue has been addressed 
by several authors in the past under different assumptions, 
and basically divides into two steps. The first involves the 
computation of the local (i.e. evaluated by an observer at 
the star surface) spectrum emitted by each patch of the star 
surface while the second requires to collect the contributions 
of surface elements which are "in view" at different rotation 
phases, making proper account for the fact that only rays 
propagating parallel to the line-of-sight (LOS) actually reach 
the distant observer. Details on each are presented in the fol- 
lowing two subsections (Sii l2.2lE3l and further in Appendix 
IXI here we discuss some general assumptions which are at 
the basis of our model. 

We take the neutron star to be spherical (mass M, ra- 
dius R) and rotating with constant angular velocity uj — 
2n/P, where P is the period. Since XDINs are slow rota- 
tors (P ~ 10 s), we can describe the space-time outside 
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the NS in terms of the Schwa rzschild geometry (see e.g. 
ICadeau. Leahy fc Morsin3l2005l for a more complete discus- 
sion about the effects of rotation). 

The star magnetic field is assumed to possess a core- 
centered dipole-|-quadrupole topology, B = Bdip + Jiquad- 
Introducing a polar coordinate system whose axis coincides 
with the dipole axis, the (polar) components of the dipole 
field at the star surface are 



Bdip,r — fdipBpCOsO 

Bdip,0 = gdipBpam9/2 



B. 



dip,4) 



= 0, 



(1) 

(2) 
(3) 



where Bp is the field strength at the magnetic pole, 9 and cj> 
are the magnetic colatitude and azimuth. The functions fdip 
and Qdip account for the effect of gravity and depend on the 
dimensionless star radius x = R/Rs wi th Rs = 2GM/c?\ 
t heir ex plicit expressions can be found in lPaee fc Sarmientd 
( 1 19961 . see also references therein) . The quadrupolar field 
can be easily derived from the spherical harmonics expansion 
and its expression, again at r = i?, can be cast in the form 



=0 



quad 



(4) 



where the qis are arbitrary constants. The polar compo- 
nents of the five g enerating vectors B^'J^^ are reported in 
IPaee &; Sarmientd il99d) . We just note that their expres- 
sion for the radial component of the zeroth vector contains 
a misprint and should read B^f^ad r ~ (3cos^ B — l)/2. Gen- 
eral relativistic effects are included in the quadrupolar field 
by multiplying the radial and angular components by the 
two functions fouad(x) a nd gquad(x) respectively (see again 
IPaee fc Sarmientoll996l for their expressions and further de- 
tails) . 

The NS surface temperature distribution, Ts, will in 
general depend on how heat is transported through the 
star envelope. Under the assumption that the field does 
not change much over the envelope scale-height, heat trans- 
port can be treated (locally) as one-dimensional. The sur- 
face temperature then depends only on the angle between 
the field and the radial directi on, cos a = B • n, and on 
the local field strength i? (see IPagelll99s|) . As shown by 
iGreenstein fc Hartkel l)l983h . an useful approximation is to 
write 



Ts — Tp 



El 



1/4 



(5) 



where the ratio of the conductivities perpendicular {K±) and 
parallel {K\\ )to the field is assumed to be constant. The polar 
value Tp fixes the absolute scale of the temperature and is 
a model parameter iPagd|l9 95 and references therein). For 
field strengths 3> 10^^ G, the conductivity ratio is much less 
than unity and eq. @ simplifies to 



1/2 



(6) 



This expression is used in the present investigation. An 
example of the thermal surface distribution induced by a 
quadrupolar field is shown in figure Q Different approaches 
which account for the variation of the conductivities (e.g. 
iHevl fc Hernauistll998ll yield similar results. Quite recently 



lGepper~ Kiiker fc Pag3 (|2004) investigated the influence of 
different magnetic field configurations on the surface tem- 
perature distribution. They find that, contrary to star- 
centered core fields, crustal fields may produce steeper sur- 
face temperature gradients. The inclusion of this effect is 
beyond the purpose of this first paper. However, we caveat 
that, for temperatures expected in XDINs (~ 10^ K), the 
differences between the two magnetic configurations start to 
be significant when the field strength is > 10^^ G. 

2.2 Radiative Transfer 

The properties of the radiation spectrum emitted by highly 
magnetized, coolin g NSs have been thoroughly discussed in 
the literature (e . g. IShibanoy et al]ll992l: IPavlov et alJll994l : 
IZane et al ."200l': Ho fc Lai"200lt lHo fc LaJl2003ft . Since the 
pressure scale-height is much smaller than the star radius, 
model atmospheres are usually computed in the plane par- 
allel approximation. Besides surface gravity, the spectrum 
emerging from each plane parallel slab depends depends 
both on the surface temperature Ts and magnetic field B, 
either its strength and orientation with respect to the local 
normal, which coincides with the unit vector in the radial 
direction n. In order to proceed we introduce a {6, (jj) mesh 
which naturally divides the star surface into a given number 
of patches. Once the magnetic field has been specified, each 
surface element is characterized by a precise value of _B, a 
and Ts- The atmospheric structure and radiative transfer 
can be then computed locally by approximating each atmo- 
spheric patch with a plane parallel slab, infinitely extended 
in the transverse direction and emitting a total flux aT*. 

Radiative tran sfer is s olved numerically using the ap- 
proach described in iLlovdl (|20o3). The code is based on the 
normal mode approximation for the radiation field propa- 
gating in a strongly magnetized medium and incorporates 
all relevant radiative processes. The full angle and energy 
dependence in both the plasma opacities and the radiation 
intensity is retained. In this respect we note that for an 
oblique field [i.e. (B/_B) • n 7^ 1] the intensity is not sym- 
metric around n and depends explicitly on both propaga- 
tion angles. If k is a unit vector along the photon direction, 
at depth r in the atmosphere the intensity has the form 
/ = Ie{t, jJ., if) where E is the photon energy, /i = n ■ k and 
ip is the associated azimuth. Calculatio ns are restricte d to 
a completely ionized H atm osphere (see iHo et aLll2003l and 
iPotekhin fc ChabriedlioO^ for a treatment of ionization in 
H atmospheres). 

Since the numerical evaluation of model atmospheres is 
computationally quite demanding, especially for relatively 
high, oblique fields and Ts < 10^ K, we preferred to cre- 
ate an archive beforehand, by computing models for pre- 
assigned values of cos a, B and Ts- The range of the lat- 
ter two parameters should be wide enough to cover the 
surface variation of B and T in all the cases of interest: 
12 < logB < 13.5 and 5.4 logT^ 6.6. Particular 
care was taken to generate models in the entire a domain, 
^5 cosQ 1. According to the adopted surface temper- 
ature distribution (eq. [S|), the regions close to the equa- 
tor have very low temperatures and can not be associated 
with any model in the archive. However, being so cool, 
their contribution to the observed spectrum is negligible 
(see >I2.4I| . For each model the emerging radiation inten- 
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sity, Ie{t = 0,/i, ((j), is stored, for 0.01 keV < -E < lOkeV, 
^ /X 1 and Q ^ if ^ 2tt. The archive itself consists of 
a six-dimensional array I{B ,Ta, cos a; E, /j., Lp) which asso- 
ciates at each set of the parameters B, cos a, Ts the (dis- 
crete) values of the angle- and energy-dependent intensity. 
Actually, since the code makes use of adaptive computa- 
tional meshes, the emerging intensities have been interpo- 
lated on common energy and angular grids before storage. 
The final array contains the emerging intensity for about 
100 atmospheric models, evaluated at 100 energy bins and 
on a 20 X 20 angular mesh, {fi, tp). 



2.3 The Observed Spectrum 

The problem of computing the pulse profile produced by hot 
caps onto the surface of a slowl y rotating NS including grav- 
ity eff ects was first tackled bv iPechenick. Ftaclas fc CohenI 
il983h . Their results were then generalized to the case 
of emission from the entire star s urface with an assigned 
tempe rature distribution by iPaed Eggs) and iLlovd et alJ 
l|2004) , in case of isotropic and non- isotropic radiation fields 
respectively. The approach used here follows that discussed 
in the two papers quoted above which we refer to for more 
details. For the sake of clarity, we will present a Newto- 
nian derivation first. Relativistic ray-bending can be then 
accounted for quite straightforwardly. 

The NS viewing geometry is described in terms of two 
angles x a^nd ^ which give the inclination of the LOS and 
of the dipole axis with respect to the star spin axis. Let z, 
hdip and p denote the unit vectors along the same three 
directions. Let us moreover introduce two cartesian coordi- 
nate systems, both with origin at the star center: the first, 
(X, Y, Z), is fixed and such that the Z-axis coincides with 
the LOS while the X-axis is in the (z, p) plane; the second, 
{x, y, z), rotates with the star. The z-axis is parallel to h^ip 
while the choice of the a;-axis will be made shortly. Each 
cartesian frame has an associated polar coordinate system, 
with polar axes along Z and z, respectively. The colatitude 
and azimuth are (Q, $) in the fixed frame, and {9, (f)) in the 
rotating one (the latter are just the magnetic colatitude and 
azimuth introduced in S I2.H . 

In the following we shall express vectors through their 
components: these are always the cartesian components re- 
ferred to the fixed frame, unless otherwise explicitly stated. 
The same components are used to evaluate both scalar 
and vector products. Upon introducing the phase angle 
7 = (jt, it follows from elementary geometrical considera- 
tions that p = (— sin x, 0, cos x) and hdip ~ {— sin x cos ^ — 
cos X sin ^ cos 7, — sin ^ sin 7, cos x cos ^ -|- sin x sin ^ cos 7) . It 
can be easily verified that q — (cos ^ sin 7, cos 7, sin ^ cos 7) 
represents an unit vector orthogonal to p and rotating with 
angular velocity u>. We then choose the a;- axis in the direc- 
tion of the (normalized) vector component of q perpendic- 
ular to hdip, 



q - {hdip ■ q) hd 



[1 - (hd^p ■ qf] 



1/2 ' 



(7) 



the j/-axis is parallel to hdip x qx. The local unit normal 
relative to a point on the star surface of coordinates (O, "!>) 
is readily expressed as n = (sin O cos "1>, sin O sin <1>, cos O). 



By introducing the unit vector nj_ , defined in strict analogy 
with q_L (see eq. [7|), the two expressions 



cos 6 
cos (h 



hdip ■ n 
n± ■ q± 



(8) 
(9) 



provide the relations between the two pairs of polar angles, 
the geometrical angles ^, x and the phase. While direct in- 
version of (|HJ poses no problems since it is $J S ^ tt, care 
must be taken to ensure that 0, as obtained from JUJ, covers 
the entire range [0, 2n]. This is achieved by replacing (j) with 
2-K — 4> when n • {hdip x qx) < 0. 

We are now in the position to compute the total 
monochromatic fiux emitted by the star and received by 
a distant observer. This is done by integrating the specific 
intensity over t he visible part of t he star surface at any given 
phase (see e.g. lLlovd et al.ll20o3l 



Fb(7) oc 



d<l> 



T{B, Tj,, cos a; E, /i, ip) du 



(10) 



where u — sinO. Further integration of eq. 1101 over 7 pro- 
vides the phase-averaged spectrum. As discussed in § 12.21 
the intensity depends on the properties of the surface patch 
and on the photon direction. The magnetic field strength 
B can be directly computed from the polar components of 
B (see § 12.11 and Page fc S armiento 1996) . The magnetic 
tilt angle a and the surface temperature (see eq. [B|) follow 
from cosa = n • B/B = Br/B, being n the unit radial vec- 
tor. The local values of B and Ts depend on {9, (j)). They 
can be then easily expressed in terms of (O, $) for any given 
phase using eqs. l|HJ-lO. Because the star appears point-like, 
there is a unique ray which reaches the observer from any 
given surface element and this implies that also fi and <p 
are functions of (O, $). It is clearly fj, = cosO, while ip is 
given by cos ip — m • v. The two unit vectors which enter 
the previous expression are, respectively, the projections of 
B and z on the plane locally tangent to the surface. They 
are expressed as m = (cos O cos cos O sin — sin O) and 
V = (B/B — n cos a) / sin a. The cartesian components of B 
needed to evaluate cos (p are derived in the Appendix. 

Gravity effects (i.e. relativistic ray-bending) can be now 
included in a very simple way. The local value of the colati- 
tude O is, in fact, related to that measured by an observer 
at infinity by the "ray-tracing" integral 

nl/2 



e 



1 - 



- 1- 



2v 



-1/2 



dv 



(11) 



where x = R/Rs- Since we are collecting the contributions 
of all surface elements seen by a distant observer, each patch 
is labelled by the two angles and $. This means that the 
integrand in IIOI I is to be evaluated precisely at the same two 
angles and is tantamount to replace O with O in all previous 
expressions. Note, however, that the innermost integral in 
noi l is always computed over Q. 

Effects of radiative beaming are illustrated in fig.|5| were 
we compare phase average spectra and lightcurves computed 
by using radiative atmosphere model with those obtained 
under the assumption of isotropic blackbody emission. As 
we can see, the pulse profiles are substantially different in 
the two cases. Also, by accounting for radiative beaming 
allow to reach relatively large pulse fractions (~ 20%). 
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2.4 Numerical Implementation 

The numerical evaluation of the phase-dependent spectrum 
has been carried out using an IDL script. Since our im- 
mediate goal is to check if the observed pulse profiles can 
be reproduced by our model, we start by computing the 
lightcurve in a given energy band, accounting for interstel- 
lar absorption and the detector response function. Most of 
the recent observations of X-ray emitting INSs have been 
obtained by XMM-Newton, so here we refer to EPIC-pn re- 
sponse function. Both absorption and the detect or respo nse 
depend upon the arrival photon energy E = Ey^l — 1/x, so 
the pulse profile in the [Ei, E2] range is given by 

Fin) oc / d$ / duM A{E)exp[-NH(T{E)] 

Jo Jo J El 

xl {B,Ts, COS a; E,fi, If) dE (12) 

f 27r pi 

d$ / du^J 
'0 Jo 

where A is the response function, Nh is the column den- 
sity, and a the interstellar a bsorption cross section (e.g. 
iMorrison fc McCammonlll98.t l. 

Since the energy integral J does not involve geometry, it 
is evaluated first. We select three energy bands, correspond- 
ing to the soft (0.1-0.5 keV) and hard (0.5-1 keV) X-ray 
colors, and to the total range 0.1-5 keV. Results are stored 
as described in S 12.21 in the case of the quantity X, the only 
difference being that energy index in the array J runs now 
from 1 to 3 in correspondence to the three energy intervals 
defined above. We then introduce a. {jj, — cosO, $) mesh by 
specifying 50 x 50 equally spaced points in the [0, 1] and 
[0, 27r] intervals, and interpolate the intensity array at the 
required values of /i, /i = u. Next, the values of O relative to 
the u grid are computed from eq. lllll . All these steps can 
be performed in advance and once for all, because they do 
not depend on the viewing geometry or the magnetic field. 

Then, once the angles Xi £, S'lid the phase 7 have been 
specified, the magnetic colatitude and azimuth, 0(0, $,7), 
<^(0, 7) can be evaluated. We use 32 phase bins and 
a set of five values for each of the two angles Xi C = 
(0° , 30° , 50° , 70° , 90° ) . The magnetic field is assigned by pre- 
scribing the strength of the quadrupolar components rela- 
tive to polar dipole field, bi — Qi/Bp {i = 1, . . . , 5), in ad- 
dition to Bp itself; in our grid, each bi can take the values 
(0, ±0.25, ±0.5). For each pair 9, (j) we then compute B and 
cos (fi; cos a gives the surface temperature Ts once Tp has 
been chosen. The corresponding values of the intensity are 
obtained by linear interpolation of the array J. Surface el- 
ements emitting a flux two order of magnitudes lower than 
that of the polar region (o-T*) were assumed not to give 
any contribution to the observed spectrum. Finally, direct 
numerical evaluation of the two angular integrals in 11211 
gives the lightcurves. Although in view of future application 
we computed and stored the lightcurves in the three energy 
bands mentioned above, results presented in § |21 and 0| are 
obtained using always the total (0.1-5 keV) energy band. 
To summarize, each model lightcurve depends on x, ^^'^ 
the five bi. No attempt has been made here to vary also Tp 
and Bp, which have been fixed at 140 eV and 6 x 10^^ G 
respectively. A total of 78125 models have been computed 
and stored. Their analysis is discussed in §|3 



3 ANALYZING LIGHTCURVES AS A 
POPULATION 

As discussed in § |^ under our assumptions the computed 
lightcurve is a multidimensional function which depends in 
a complex way on several parameters. Therefore, an obvi- 
ous question is whether or not we can identify some possible 
combinations of the independent parameters that are asso- 
ciated to particular features observed in the pulse shape. 
The problem to quantify the degree of variance of a sam- 
ple of individuals (in our case the lightcurves) and to iden- 
tify groups of "similar" individuals within a population is 
a common one in behavioral and social sciences. Several 
techniques have been extensi vely detailed in many books 
of muhivariate statistics (e.g. lKendaIilll957l: iManlvl Il998h 
and, although they have been little used in physical sci- 
ences, a few applications to astrophysical problems have 
been presented over the past decades (see e.g.'Whitnev'l983; 
iMittaz. Penston fc Sniiderj [r990: Hcver & Schloerb 1997*). 

We focus here on a particular tool called principal 
components analysis (PCA), which appears promising for 
a quantitative classification of the lightcurve features. The 
main goal of PCA is to reduce the number of variables that 
need to be considered in order to describe the data set, by 
introducing a new set of variables Zp (called the principal 
components, PCs) which can discriminate most effectively 
among the individuals in the sample, i.e. the lightcurves in 
our present case. The PCs are uncorrelated and mutually- 
orthogonal linear combinations of the original variables. Be- 
sides, the indices of the PCs are ordered so that zi dis- 
plays the largest amount of variation, Z2 displays the second 
largest and so on, that is, var(2i) ^ var(z2) ^ . • . ^ var(zp) 
where var(2:fe) is the variance in the sample associated with 
the fc-th PC. Although the physical meaning of the new 
variables may be in general not immediate, it often turns 
out that a good representation of the population is obtained 
by using a limited number of PCs, which allows to treat 
the data in terms of their true dimensionality. Beyond that, 
PCA represents a first step towards other kind of multivari- 
ate analyses, like the cluster analysts. This is a tool which 
allows to identify subgroups of objects so that "similar" ones 
belong to the same group. When applying the cluster anal- 
ysis algorithm a PCA is performed first in order to reduce 
the number of original variables down to a smaller number of 
PCs. This can substantially reduce the computational time. 

Since both tools are extensively described in the liter- 
ature, we omit all mathematica l details f or wh ich a n inter- 
ested reader may refer to, e.g.. iKendalJ lll957ll and iManTvl 
il998h . Let us denote with j/ps (p = 1, . . . , P; s = 1, . . . , 5) 
the values of the intensity computed at each phase bin, for 
each model lightcurve. Let us also introduce the "centered" 
data Xps, defined as 



— iUpS /ip) /Sp , 



(13) 



where /ip and Sp are the mean and the variance of the com- 
puted data, respectively. In order to shift to the new PCs 
variables Zpa, we computed the transformation matrix V' 
such that 



El 
VpqXqs ■ 



(14) 



A sufficient condition to specify univocally V' is to impose 
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that the axes of the new coordinate system (i.e. the PCs) 
are mutually orthogonal and linearly independent. 

By applying the PCA to our sample of models, we have 
found that each lightcurve can be reproduced using only the 
first 20-21 more significant PCs (instead of 32 phases) and 
that the first 4 PCs alone account for as much as ~ 85% 
of the total variance. Moreover, ~ 72% of the variance is 
in the first three PCs only. It is therefore meaningful to 
introduce a graphical representation of the model dataset in 
terms of the first three Zi's. This is shown in figure |3 where 
black/red squares gives the position in the ziZ2Zz space of 
quadrupolar/dipolar models. To better visualize the latter, 
an additional set of lightcurves was computed, bringing the 
total number of dipolar models displayed in fig.j^to 100. 

An insight on the lightcurve property each of the PCs 
measures can be obtained by inspecting the coefficients 
v'pq of the linear combination which gives the Zp's for an 
assigned dataset [see eq. I|14|l ]. Since Zp — '^^'VpqXq oc 

Jq^ Up(7)-F(7) d'-f, this is tantamount to assess the effect of 
the filter Vp{'y) on the lightcurve F{'y) knowing the values 
of the former at a discrete set of phases. The first four Vp 
are shown in Fig. |1] The first PC provides a measure of the 
amplitude of the pulse; it is always zi > 0, and large val- 
ues of zi correspond to low pulsed fractions. Both Z2 and 
Z3 may take either sign (same for higher order PCs) and 
give information on the pulse shape. We note that, although 
the absolute phase is in principle an arbitrary quantity, the 
whole models population has been computed using the same 
value. Therefore, when studying the morphological proper- 
ties of the sample of lightcurves, it is meaningful to refer to 
the symmetry properties with respect to the parameter 7. 
Large and negative values of Z2 imply that the main con- 
tributions to the lightcurve comes from phases around zero. 
Z3 measures the parity of the lightcurve with respect to half 
period. Pulses which are symmetric have 2:3 — 0. As fig. 
131 clearly shows, the PCA is very effective in discriminating 
purely dipolar from quadrupolar models. The former cluster 
in the "tip of the cloud" , at large values of Zi, negative values 
of Z2 in the 23 = plane, as expected. In fact, dipolar pulse 
patterns are always symmetrical and their pulsed fraction 
is quite small (semi-amplitude < 10%). It is worth noticing 
that quadrupolar magnetic configurations too can produce 
quite symmetrical lightcurves, e.g. the black squares in fig. 
I^with Z3 ~ 0. However, in this case the pulsed fraction may 
be much larger, as indicated by the lower values of zi at- 
tained by quadrupolar models with 23 = in comparison 
with purely dipolar ones. This implies that a symmetrical 
lightcurve is not per se an indicator of a dipolar magnetic 
geometry. 

As suggested by some authors (e.g. lHeckll976l : fWhitnevl 

^S^), PCs can then be used as new variables to describe 
the original data. However, in the case at hand, the prob- 
lem is that although PCs effectively distinguish among pulse 
patterns, they have a purely geometrical meaning and can 
not be directly related to the physical parameters of the 
model (&i,XiC)- We have found that the standard regres- 
sion method does not allow to link the PCs with the model 
parameters, which is likely to be a signature of a strong non- 
linearity (see §|5]|. Instead, the PCA can be regarded as a 
method to provide a link between the numerous lightcurves, 
in the sense that models "close" to each other in the PCs 



space will have similar characteristics. Unfortunately, de- 
spite different definitions of metrics have been attempted, 
so far we found it difficult to translate the concept of "prox- 
imity" in the PCs space in a corresponding "proximity" in 
the 7-dimensional space of the physical parameters 5, x a-nd 
bi {i — 1, ... ,5). By performing a cluster analysis we found 
that two separate subgroups are clearly evident in the PCs 
space, one of which encompasses the region occupied by 
purely dipolar models (see fig. El top panel). However, again 
it is not immediate to find a corresponding subdivision in 
the physical space. Due to these difficulties, we postpone a 
more detailed statistical analysis to a follow-up paper, and, 
as discussed in the next section, we concentrate on the di- 
rect application of our models to the observed lightcurves of 
some isolated neutron stars. 



4 AN APPLICATION TO XDINSS 

In the light of the discussion at the end of the previous sec- 
tion, the PCA may be used to get insights on the ability of 
the present model to reproduce observed lightcurves. A sim- 
ple check consists in deriving the PC representation of the 
pulse profiles of a number of sources and verify if the cor- 
responding points in the 2:12:223 space fall inside the volume 
covered by the models. We stress once again that, although 
the model lightcurves depend upon all the PCs, zi, 22, and 
2:3 alone provide a quite accurate description of the dataset 
since they account for a large fraction of the variance. In 
this sense the plots in fig. 01 give a faithful representation of 
the lightcurves, with no substantial loss of information: pro- 
files close to each other in this 3-D space exhibit a similar 
shape. To this end, we took the published data for the first 
four pulsating XDINSs listed in table and rebinned the 
lightcurves to the same 32 phase intervals used for the mod- 
els. Since the PCA of the model dataset already provided 
the matrix of coefficients v'pq (see eq. 1141 '). the PCs of any 
given observed lightcurve are Zp'"' — ^ v'pqXq'"' , where Xq^" 
is the (standardized) X-ray intensity at phase 7,. As it can 
be seen from fig. the observed pulse profiles never fall in 
the much smaller region occupied by purely dipolar models. 
However, they all lie close to the quadrupolar models, indi- 
cating that a quadrupolar configuration able to reproduce 
the observed features may exist. A possible exception to this 
is the first observation of RX J0720.4-3125 {XMM-Newton 
rev. 078), for which the pulse profile appears quite symmet- 
ric and the pulsed fraction is relatively small (see table 0. 
While a purely dipolar configuration may be able to repro- 
duce the lightcurve for rev. 078, a visual inspection of fig.jHl 
shows that this is not the case for the second observation of 
the same source {XMM-Newton rev. 711). Despite the pulse 
shape is still quite symmetrical, the pulsed fraction is defi- 
nitely too large to be explained, within the current model, 
with a simple dipole field. The same considerations apply to 
the lightcurve of RX J0420. 0-5022. 

Just as a counter-example, we added in fig. (21 the PC 
representation of the X-ray lightcurve of the Anomalous X- 
ray pulsar IE 1048.1-5937 observed in two differe nt epochs, 
2000 December and 2003 June fsee lMereehetti e t al. 2001). 
The new data points fall well outside the region of the popu- 
lated by quadrupolar models and can be seen only in the last 
panel of fig. (21 In the case of this source, the pulsed fraction 
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is so high (89% and 53% in June and December respectively) 
that no quadrupolar configuration could account for it. In 
terms of PCs, both points have a relatively low value of zi 
{zi = 9.9 and zi = 6.9). In this case no fit can be found, not 
surprisingly, since we do expect a large contribution from a 
non-thermal power law component to the X-ray emission of 
Anomalous X-ray pulsars. 

To better understand to which extent quadrupolar sur- 
face distributions may indeed produce lightcurves close to 
the observed ones, we select the model in the data set which 
is "closer" to each of the observed lightcurves. This is done 
looking first for the minimum Euclidean distance between 
the observed and the model pulse profiles in the PCs space. 
Note that in this case all the relevant PCs were used, that 
is to say the distance is defined by = ~ Zp*"^)^. 

The computed model which minimizes D is then used as 
the starting point for a fit to the observed lightcurve which 
yields the best estimate of the model parameters. The fitting 
is performed by computing "on the fiy" the lightcurve for 
the current values of the parameters and using the standard 
(i.e. not the PC) representation of both model and data. The 
quadrupolar components and viewing angles are treated as 
free parameters while the polar values Tp and Bdip are fixed 
and must fall in the domain spanned by our atmosphere 
model archive.^ The (arbitrary) initial phase of the model 
is an additional parameter of the fit. 

The results of the fits are shown in figures |^ |S| and [7| 
for Bdip = 6 X 10^^ G, logrp(K) = 6.1 - 6.2. It is apparent 
that the broad characteristics of all the XDINSs lightcurves 
observed so far may be successfully reproduced for a suit- 
able combination of quadrupolar magnetic field components 
and viewing angles. However, although in all cases a fit ex- 
ists, we find that in general it is not necessary unique. This 
means that the model has no "predictive" power in provid- 
ing the exact values of the magnetic field components and 
viewing angles. For this reason we do not attempt a more 
complete fit, i.e. leaving also Tp and Bdip free to vary, nor 
we derive parameter uncertainties or confidence levels. Our 
goal at this stage has been to show that there exist at least 
one (and more probably more) combination(s) of the param- 
eters that can explain the observed pulse shapes, while this 
is not possible assuming a pure dipole configuration. 

The case of RX J0720. 4-3125 deserves, however, some 
further discussion. This source, which was believed to be 
stationary and as such included among XMM-Newton EPIC 
and RGS calibration sources, was recently found to exhibit 
rather sens ible variations both in its spectral an d timing 
properties llDe Vries et al.1 120041: fvink et al.ll2004) . In par- 
ticular, the pulse shape changed ove r the 2 yrs time be - 
tween XMM revolutions 78 and 711. iDe Vries et al.l ll2004ll 
proposed that the evolution of RX J0720. 4-3125 may be pro- 
duced by a (freely) precessing neutron star. This scenario 
can be tested by our analysis, since in a precessing NS only 
the two angles ^ and x a-re expected to vary while the mag- 
netic field remains fixed. This means that having found one 
combination of parameters which fits the lightcurve of rev. 



^ In general this will not contains the exact values inferred from 
spectral observations of XDINSs. However, we have checked that a 
fit (albeit with different values of the quadrupolar field) is possible 
for different combinations of Bdip a-nd Tp in the range of interest. 



78, a satisfactory fit for rev. 711 should be obtained for the 
same 6i's and different values of the angles. Despite several 
attempts, in which the proper values of Tp as derived from 
the spectral analysis of the two observations were used, we 
were unable to obtain a good fit varying the angles only 
(see figure ISJ. We performed also a simultaneous fit to both 
lightcurves starting from a general trial solution and asking 
that the bi's are the same (but not necessarily coincident 
with those that best fit data from rev. 78), while the two 
angles (and the initial phases) are kept distinct (see figure 
1^. Both approaches clearly indicate that a change in both 
sets of quantities (magnetic field and angles) is required (as 
in Fig. [7J . From a physical point of view it is not clear how 
magnetic field variations on such a short timescale may be 
produced, therefore at present no definite conclusion can be 
drawn. For instance, one possibility that makes conceivable 
a change of the magnetic field structure and strength on a 
timescale of years is that the sur face field is small scaled 
JCeppert. Rheinhardt fc"Giill2003l) . In this case, even small 
changes in the inclination between line of sight and local 
magnetic field axis may cause significant differences in the 
"observed" field strength. 



5 DISCUSSION 

X-ray dim isolated neutron stars (XDINs) may indeed rep- 
resent the Rosetta stone for understanding many physical 
properties of neutron stars at large, including their equation 
of state. Despite their potential importance, only recently 
detailed observations of these sources have become progres- 
sively available with the advent of Chandra and XMM- 
Newton satellites. These new data, while confirming the 
thermal, blackbody-like emission from the cooling star sur- 
face, have revealed a number of spectral and timing features 
which opened a new window on the study of these objects. 
Some issues are of particular importance in this respect: i) 
the discovery of broad absorption features at few hundreds 
eVs, ii) the quite large pulsed fractions, iii) the departure 
of the pulse shape from a sine wave, and iv) the long-term 
evolution of both the spectral and timing properties seen in 
RX J0720.4-3125 and, to some extent, in RX J0420.0-5022. 
Pulse-phase spectroscopy confirms that spectral and timing 
properties are interwoven in a way which appears more com- 
plex than that expected from a simple misaligned rotating 
dipole, as the anti-correlation of the absorption line strength 
and of the hardness ratio with the intensity along the pulse 
testify. 

Motivated by this, we have undertaken a project aimed 
at studying systematically the coupled effects of: i) radia- 
tive transfer in a magnetized atmospheric layer, ii) thermal 
surface gradients, and iii) different topologies of the surface 
magnetic field in shaping the spectrum and pulse shape. 
The ultimate goal of our investigation is to obtain a simul- 
taneous fit of both pulse profile and (phase-averaged and 
-resolved) spectral distribution. As detailed comparisons of 
synthetic spectra with observations have shown, no com- 
pletely satisfactory treatment of spectral modelling for these 
sources is available as yet. For this reason, here we presented 
the general method and concentrated on the study of the 
lightcurves, computed assuming a pure H, magnetized at- 
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mosphere . The pulse shapes, in fact, should be less sensitive 
on the details of the chosen atmospheric model. 

We caveat the reader that our results have been com- 
puted under a number of simplifying assumptions. For in- 
stance, there are still considerable uncertainties in current 
modelling of the NS surface thermal emission: we just men- 
tion here that most of the observed NS spectra cannot be 
reproduced by the theoretical models currently available 
(see lHaber]||200i iHaberll l20o3 for reviews and references 
therein). Second, the surface temperature has been derived 
using eq. lO which is based on the assumption that the 
temperature profile is only dictated by the heat transferred 
from the star interior. While this is expected to be the main 
mechanism, other effects may significantly contribute. For 
instance, heating of the star surface may occur because of 
magnetic field decay, and the polar caps may be re-heated 
by back-flowing particles or by internal friction. Third, we 
have computed the emissivity by assuming that each atmo- 
spheric, plane parallel, patch emits a total flux aT^. In other 
words, while the spectral distribution is computed using a 
proper radiative transfer calculation, we introduced the fur- 
ther assumption that each slab emits the same total flux as 
a blackbody radiator. A more consistent approach would be 
to compute the spectrum emitted by each patch by taking 
the value of Ts and the efficie ncy of the crust as a radiator 
jTurolla. Zane fc Drakd ''2004') as the boundary condition at 
the basis of the atmosphere. Our working assumption avoids 
the burden of creating a complete grid of models in this pa- 
rameter, with the disadvantage that the spectral properties 
of each patch may be slightly approximated. 

As far as the application to XDINSs is concerned, the 
greatest uncertainties arise because in this paper we are not 
fitting simultaneously the observed spectra and pulse shape. 
The quadrupolar components and viewing angles are treated 
as free parameters while the polar values Tp and B^ip are 
fixed and, of course, they must fall in the domain spanned by 
our atmosphere model archive. As stated earlier, albeit we 
have checked that a fit is possible for different combinations 
of Bdip and Tp in the allowed range, in general the archive 
will not contain the exact values of T inferred from spectral 
observations of the coldest XDINSs. As long as we make the 
same assumption that the local spectrum emitted by each 
patch of the star surface is computed using fully ionized, 
magnetized hydrogen atmosphere models, we do still expect 
to reach a good fit by using a temperature of a few tens 
of eV smaller than those used, albeit with different values 
of the quadrupolar components. However, partial ionization 
effects are not included in our work, and bound atoms and 
molecules can affect the results, changi ng the radiation prop- 
erties at relatively low T a nd high B jPotekhin fc Chabrierl 
l2003l:lPotekhin et aljliooi) . 

Even more crucial is the issue of fitting with a realistic 
value of the magnetic field. The recent detection of absorp- 
tion features at « 300-700 eV in the spectrum of XDINSs 
and their interpretation in terms of proton cyclotron reso- 
nances or bound-bound transitions in H or H-like He, may 



indicate that these sources possess strong 


; magnetic fields, up 


to B - 9 X 10^^ G fVan Kerkwiik et al. 


^2004; Habcrl ctM 



l2004l:IZane et al.ll2005l) .' A (slightly) more direct measure- 
ment, based on the spin-down measure, has been ob- 
t ain ed so far only in one source (i.e. RX J0720.4-3 125, 
e.g. ICropper et aLll2004 iKapl an fc van Kerkwiiklliooi) . In 



this case the measurement points at a more moderate field, 
B ~ (2 — 3) X 10^'' G, which is only a few times larger 
than that used here. However, the possibility that (some) 
of these objects are ultra-magnetized NSs is real. Would 
this be confirmed, our model archive should be extended to 
include higher field values. This, however, poses a serious 
difficulty, since the numerical convergence of model atmo- 
spheres is particularly problematic at such large values of 
B for T < lO'' K. Moreover, as mentioned in t|2.1l if such 
high field strengths are associated to crustal (and not to 
star-centered) fields, the surface temperature gradient is ex- 
pected to be substan tially different from that used in the 
present investigation iGeppert. Kiiker fc Pagell2004ll . 

The application presented here makes only use of the 
properties of the pulse profile in the total energy band, and 
does not exploit color and spectral information available for 
these sources. This worsens the problem of finding a unique 
representation, problem which is to some extent intrinsic, 
due to the multidimensional dependence of the fitting func- 
tion on the various physical parameters. Aim of our future 
work is to reduce the degeneracy by refining the best fit solu- 
tions by using information from the light curves observed in 
different color bands and/or from the spectral line variations 
with spin pulse. Also, a more detailed statistical analysis on 
the models population based on algorithms mor e sophisti- 
cated than the PCA's (lGifilll99d : ISaegusal l2005') may shed 
light on the possible correlation between the physical param- 
eters in presence of strong non-linearity and possibly on the 
meaning of the two subclasses identified through a cluster 
analysis. 
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Isolated Neutron Stars Parameters. 



Source P (s) Semi-Ampl. XMM-Newton rev. Refs. 



RX J0420.0-5022 3.45 13% 570 1 

RX J0806. 4-4123 11.37 6% 618 1 

RBS 1223 10.31 18% 561 2 

RX J0720.4-3125 8.39 11% 78 3 

RX J0720.4-3125 8.39 16% 711 3 

RBS 1774 9.44 4% 820 4 



Table 1. The light curves used in this paper are taken fr om references listed in the fourth column, which correspond to: fll lHaberl et alj 
Jiooil); f2l lHaberl et al.U2nO.'J^ : [3l lDe Vries et al.l J2004h : l4l lZane et aIli20n,'J) . 




Figure 1. Hammer projection illustrating tlie thermal surface distribution of a NS; the quadrupolar field components (in units of the 
polar dipolar field) are 60 = 0.0248, 61 = 0.684, 62 = 0.0501, 63 = 0.0480, 64 = -0.0624. The contours are labelled by the values of 
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Figure 2. Effects of Radiative beaming. In eacli pair of figures we show pliasc average spectrum (left) and light curves in the 0.1-2 kcV, 
0.1-0.5 kcV and 0.5-2 kcV bands (riglit, arbitrary normalization). In all panels solid lines refer to atmospheric models, while dashed lines 
to blackbody emission. The four pairs of figures correspond to different values of g : from top left to bottom right it is ^ = 0° , 30° , 60° , 90° , 
respectively. The remaining parameters are fixed at: B^ip = 6x10^'^ G, Tpoi = 2.5MK, bo = 0.5, 6i = 63 = 64 = 0, 62 = 0.9, x = 90°- 
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This figure is provided separately as the f3.png file 



Figure 3. The computed population of lightcurves plotted against the first 3 principal components; units on the axes are arbitrary. 
Yellow symbols mark the position of the EPIC-PN light curves of XDINSs (both observations relative to rev. 78 and rev. 711 arc shown 
for RX J0720. 4-3125), red symbols mark the lightcurves computed assuming a pure dipolar field. The PCs representation (limited to the 
first three PCs, which alone account for the 72% of the total variance) of the observed lightcurves falls within the domain spanned by the 
quadrupolar model representations; this is why at least one fitting solution can be found. For comparison, the position of the EPIC-PN 
lightcurves of the Anomalous X-ray pulsar IE 1048.1-5937 is also shown. These two points fall outside the model distribution and they 
can be seen only in the last panel (violet diamonds; see text for details). 
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Figure 5. Top: Results of a cluster analysis performed on the 
models population. As in fig. |3] axis are the first three PCs: the 
different scale is due to the fact that here the PCs have been 
computed on the centered data. Two separate subgroups, plotted 
in red and blue respectively, are evident in the PCs space. Bottom: 
Fit of the EPIC-PN (0.12 -0.7 keV) lightcurve of RX J0420.0-5022 
detected during rev. 570 iHaberl et al .12004) . Data point refer to 
the smoothed observed lightcurve; dashed line: trial solution as 
inferred from the closest model in the PCs space; solid line: best 
fit solution. The best fit parameters are: bo = —0.48, bi = 0.02, 
62 = -0.25, 63 = 0.35, 64 = -0.20, S, = 39.9°, x = 91-2°. 
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Figure 6. Same as in fig.inifor two further XDINSs. Top: the 
EPIC-PN (0.12- 1.2 keV) lightcurve of RX J0806.4-4123 detected 
during rev. 618 iHaberl et all2004^ . The best fit parameters are: 
60 = 0.39, bi = -0.37, 62 = 0.12, 63 = -0.13, bi = 0.49, g = 0.0°, 
X = 59.2°. Bottom: the EPIC-PN (0.12-0.5 keV) ligh tcurve of 
RBS 1223 detected during rev. 561 <Haberl et alJl2003^ . Best fit 
parameters are: 60 = 0.21, 61 = 0.02, 62 = 0.59, 63 = 0.53, 
64 = 0.50, ? = 0.0°, X = 95.1° 
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Figure 7. Same as in fig. HJ for two different EPIC-PN (0.12- 
1.2 keV) lightcurves detected from RXJ0720.4-3125 (botii from 
iDe Vries et aill20n4ft . Top (rev. 78): bo = 0.36, 6i = 0.43, 62 = 
-0.16, bs = -0.16, 64 = -0.39, (, = 0.0°, x = 68.1°. Bottom 
(rev. 711): 60 = 0.45, 61 = 0.49, 62 = -0.06, bs = -0.08, 64 = 
-0.26, 5 = 0.0°, X = 87.6°. The second fit is obtained by leaving 
all parameters free with respect to the solution shown in the top 
panel. 



Figure 8. An attempt to fit the EPIC-PN (0.12-1.2 keV) 
Hghtcurves detected from RXJ0720. 4-3125 during rev. 711 by us- 
ing a fitting function with either the two viewing angles (top) or 
the values of bi's (i = 0...4, bottom) kept fixed and equal to the 
values derived in the case of rev. 78. The dashed line is the same 
as in Fig. |7| lower panel. In this particular example, we also al- 
lowed a fractional change of ~ 20% in the temperature between 
the two revolutions. Both attempts give an unsatisfactory result. 
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Figure 9. A simultaneous fit to both lightcurves detected from 
RXJ0720. 4-3125 during rev. 78 (shown in the phase interval 0-1) 
and rev. 711 (phase 1-2). Only changes in the viewing angles have 
been allowed between the two epochs. The result is unsatisfactory. 
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APPENDIX A: 

Let us write the magnetic field at any given point on the 
star surface in terms of its cartesian components relative to 

the rotating frame (x, y, z) 

B = B.^c\^ + By{haip x qx) + B^hup . (Al) 

The three unit vectors which identify the cartesian axes in 
the rotating frame have polar components r, 6, (relative 
to the same frame) 

qx = (sin 6 cos 8, cos 6 cos (f>, — sin 9 sin 0) 

bdip X qx = (sin 6* sin 6*, cos 9 sin 0, sin ^ cos 0) (-^-2) 

bdip = (cos 9, — sin 9, 0) . 

Taking the scalar product of B = {Br, Bg,B^) with each of 
the three unit vectors above, we have 

Bx = Br sin 6 cos 9 + Be cos 9 cos 4> — B^ sin 9 sin cj) 

By = Br sin ^ sin ^ + Be cos ^ sin + B^ sin 9 cos 4> (A3) 

Bz = Br cos 9 — Besm.9 . 

The cartesian components of B in the fixed frame are then 
given by 

Bx = Ba:{ci±)x + By{hdzp X q±)x + -Bz(bdip)x 

By = B^{q_±)Y + By{hdip X q_±)Y + B^ (b^^p ) y ( A4) 

Bz = Bxiq±)z + By{hdip x qx)z + Bz(hdip)z ■ 
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